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Abstract 

Starting from inhomogeneous time scaling and linear decorrelation between successive price re- 
turns, Baldovin and Stella recently proposed a way to build a model describing the time evolution 
of a financial index. We first make it fully explicit by using Student distributions instead of power 
law-truncated Levy distributions; we also show that the analytic tractability of the model extends 
to the larger class of symmetric generalized hyperbolic distributions and provide a full computation 
of their multivariate characteristic functions; more generally, the stochastic processes arising in 
this framework are representable as mixtures of Wiener processes. The Baldovin and Stella model, 
while mimicking well volatility relaxation phenomena such as the Omori law, fails to reproduce 
other stylized facts such as the leverage effect or some time reversal asymmetries. We discuss how 
to modify the dynamics of this process in order to reproduce real data more accurately. 



*Electronic address: damien. challet@unifr.ch 

^Electronic address: ppeirano@libero.it, corresponding author 



1 



I. HOW SCALING AND EFFICIENCY CONSTRAINS RETURN DISTRIBUTION 



Finding a faithful stochastic model of price time series is still an open problem. Not 
only should it replicate in a unified way all the empirical statistical regularities, often called 
stylized facts, (cf e.g. Bouchaud and Potters [15], Cont [21]), but it should also be easy to 
calibrate and analytically tractable, so as to facilitate its application to derivative pricing and 
financial risk assessment. Up to now none of the proposed models has been able to meet all 
these requirements despite their variety. Attempts include ARCH family (Bollerslev et al. 
[10], Tsay [50] and references therein), stochastic volatility (Musiela and Rutkowski [41] 
and references therein), multifractal models (Bacry et al. [1], Borland et al. [13], Eisler and 
Kertesz [27], Mandelbrot et al. [39] and references therein), multi-timescale models (Borland 
and Bouchaud [12], Zumbach [54], Zumbach et al. [56]), Levy processes (Cont and Tankov 
[22] and references therein), and self-similar processes (Carr et al. [18]). 

Recently Baldovin and Stella (B-S thereafter) proposed a new way of addressing the 
question. We advise the reader to refer to the original papers Baldovin and Stella [4, 5, 6] 
for a full description of the model as we shall only give a brief account of its main underlying 
principles. Using their notation let S{t) be the value of the asset under consideration at time 
t, the logarithmic return over the interval [t, t + 6t] is given by rt^st = In S{t + 6t) — In S(t); 
the elementary time unit is a day, i.e., t = 0,1,... and 6t = 1,2,... days. In order to 
accommodate for non-stationary features, the distribution of is denoted by PtMij) which 
contains an explicit dependence on t. The most impressive achievement of B-S is to build 
the multivariate distribution PQ"^(ro,i, . . . ,''"n,i) of n consecutive daily returns starting from 
the univariate distribution of a single day provided that the following conditions hold: 

1. No trivial arbitrage: the returns are linearly independent, i.e. £'(rj i,rj i) = for 
i 7^ j, with the standard condition -E'(rj i) = 0. 

2. Possibly anomalous scaling of the return distribution with respect to the time interval 
5t, with exponent D: 



3. Identical form of the unconditional distributions of the daily returns up to a possible 
dependence of the variance on the time t, i.e. 
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As shown in the addendum of Baldovin and Stella [5] these conditions admit the solution 



f^^ih, ...,K) = ~g{sjarkl + ■■■ + al^kl), (1) 

where /g"'* is the characteristic function of -Po"\ g the characteristic function of Pq,!? 
a^^ = — (z — 1)^^. In this way the full process is entirely determined by the choice of 
the scaling exponent D and the distribution Po,i- Therefore the characteristic function of 

ftAk) = /ff (o_^,^_^, 0, . . . , 0) = ~g{k^{t + 6ty^-t^^), 

t terms St terms 



I.e. 



Rst{r) = —======P,, 



The functional form of g in Eq. (1) introduces a dependence between the unconditional 
marginal distributions of the daily returns by the means of a generalized multiplication (8) 
in the space of characteristic functions, i.e., 



fo^iki, ...,kn) = g{a^ki) ®g---®g g{a^kn) 



with ®g defined by 



X I 



y = g (/rWTFW) • (2) 



At first sight this last equation may seem a trivial identity, but it does hide a powerful 
statement. Suppose indeed that instead of starting with the probability distribution g, one 
takes a general distribution with finite variance cr^ = 2 and characteristic function pi, then 
it is shown in Baldovin and Stella [4] that 

N^ooP' ii^) ^^---^SP^ (-4^) = (3) 



terms 

This means that in this framework the return distribution at large scales is independent 
of the distribution of the returns at microscopic scales: it is completely determined by 
the correlation introduced by the multiplication with fixed point g. Note that if g is 
the characteristic function of the Gaussian distribution, then (g)^ reduces to the standard 
multiplication and one recovers the standard Central Theorem Limit. 
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As the volatility of the model shrinks in an inexorable way, Baldovin and Stella propose 
to restart the whole shrinking process after a critical time Tc long enough for the volatility 
autocorrelation to fall to the noise level. In this way one recovers a sort of stationary time 
series when their length is much greater than Tc- In this case one expects that the empirical 
distribution of the return Pst{r) over a time horizon 6t <ti evaluated with a sliding window 
satisfies 



In the original papers no market mechanism is proposed for modeling the restart of the 
process; it is simply stated that the length of different runs and the starting points of the 
processes could be stochastic variables. In their simulations the length of the processes was 
fixed to r = 500, which corresponds to slightly more than two years of daily data. 

II. A FULLY EXPLICIT THEORY WITH STUDENT DISTRIBUTIONS 

In Baldovin and Stella [5] a power law truncated Levy distribution is chosen to describe 
the returns 



In Sokolov et al. [47] it is shown that this expression is indeed the characteristic function 
of a probability density with power law tails whose exponent is exponent 5 — a. How- 
ever, this choice is problematic in two respects: its inverse Fourier cannot be computed 
explicitly, which prevents a fully explicit theory. In addition, for Eq. (1) to be consistent. 



for all n. In Baldovin and Stella [5] only numerical checks are performed to verify this prop- 
erty. But as discussed for example in Bouchaud and Potters [15] both truncated Levy and 
Student distributions yield acceptable fits of the returns on medium and small time scales. 
In the present context, the Student distribution, sometimes referred to as g-Gaussian in the 
case of non-integer degrees of freedom, is a better choice; it provides analytic tractability 
while fitting equally well real stock market prices (see alsoOsorio et al. [44]). The fit of the 
daily returns of the S&P 500 index in the period with a Student distribution 




(4) 




(5) 



9(^/^1 + ■ ■ ■ + ^n) must be the characteristic function of a multivariate probability density 



r(- + -) 1 
7rV2Ar(|) (l + f)f+^ 
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Figure 1: Centered distribution of the 14956 daily returns of the S&P 500 index (January, 
3th 1950 - June, 11th 2009), and the corresponding fitting with Student {u = 3.21, 
A = 0.0109) and Gaussian distribution (a = 0.0095). 

is reported in Fig. 1[57]. 

The characteristic function of the Student density is 

~g{k) = ^^k^K^{k), (6) 

where Ka is the modified Bessel function of third kind. As demonstrated in the appendix, the 
inverse Fourier transform of g{>y kf + + k"^) for any integer n is simply the multivariate 
Student distribution (see also Vignat and Plastino [52]). The general form of this distribution 
can be written as 

o(^)fx A) = ^^^^^^ I m 

9n ^n,/2(det A)V2r(|) (1 + X* A'^x) ^ +t ' ^ ' 

where > 1 is the exponent of the power law of the tails, V{r > R) oc l/W^ and A is a 
positive definite symmetric matrix governing the variance-covariance matrix -E(xj, Xj) = 
which does exist provided that u > 2. 

In passing, the same properties are shared by multivariate symmetric generalized hyper- 
bolic distributions introduced in finance by Eberlein and Keller [26] (see also Bingham and 
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Kiesel [8]). The general case is obtained by an affine change of variable, but for the sake of 
brevity let us restrict to 

a 2 1 



(27r) 2 A|(a) (1 + r2)4 + 4 2 2 

for X G M" and r the usual euclidean norm of x. Student distributions are recovered in the 
limit a ^ 0^. As shown in the appendix, its characteristic function is given for any n by 



-fvii(a) a2 



with k = VELi ^f- 

In the following we restrict the discussion to the Student distributions. Hence we assume 
that the distribution of the return is given by Eq. (7) with characteristic function given by 
Eq. (6), where A is a diagonal matrix 



k = Vk*Ak = A-^A;2 + (22^ - l)kf + ■■■ + {n^^ - {n - iy^)kl^^ 

and governs the variance of the returns on the time scale chosen as a reference. Thanks 
to the fact that the diagonal elements of A form a telescoping series the process is indeed 
consistent for any number of discrete steps. Moreover it can be generalized to the continuous 
time by setting, in the same consistent way, 

= 9^:\roAto^n,At.^ . . . , ri„_„At„_„ A = diag(t?^, tf - t?^, . . . , - t^^i)), (8) 

where tj = Xlto ^^i^ 3 ^ ^ now A = diag(tP, tf^ -tf,..., t^^ -tl^i). The existence 
of the continuum process is then guaranteed by the Kolmogorov extension theorem. Starting 
from this expression a wider class of processes can be generated by suitable transformations 
of the time, i.e., by substituting the function ti tf^ for any monotonically increasing 
continuous function ti T(ti). The process followed by the price x{t) = InS'(t) is a Student 
process too, with same exponent u and non diagonal matrix Aij = {—iy'^^T{tmin{i,j))- 

The Student setting makes easier to interpret the correlations induced by the pointwise 
non-standard product of (2) in the characteristic function space. If we consider two variables 
Xi and X2 distributed according to gi{x), the joint probability function will be g2{xi,X2). 
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(a) 3D perspective. (b) Level plot. 

Figure 2: Student copula density with u = 3 and trivial correlation matrix. 

The variables Xi = G{xi) = J^'^dx gi{x) are distributed uniformly on the interval [0, 1]; by 
definition, the copula function c{Xi,X2) (cf. e.g. Nelsen [43] for a general theory) is 

^(Y v^- (r^~^(v ^ r^-V^ _ g2{G-\Xi),G~\X2)) 

c{X,,X2)-g2{G {X,),G ™)^^ - g^G-^^x,)) g{G-^{X2)y 

In our case c is none other than the Student copula function, generally applied in finance for 
describing the correlation among asset prices (Cherubini et al. [20], Malevergne and Sornette 
[38]). A picture of this copula density with u = 3 and A the identity matrix is given in 
Fig. 2. Although Student and generalized hyperbolic distributions are usually adopted for 
modeling returns of several assets over the same time intervals, the framework proposed by 
Baldovin and Stella allow them to model the returns of a single asset over different time 
intervals. 
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III. THE BALDOVIN-STELLA PROCESS AS MULTIVARIATE NORMAL VARI- 
ANCE MIXTURES 



According to the B-S framework we have to look for functions </> : M ^ C, such that 
(jn '■ M" — C with gn{ki, k2, . . . , kn) = (t>{k\ + /cf + ■ ■ ■ + /c^) is the characteristic function of a 
probability distribution for any n. Then from Eq. (8) we obtain a unique stochastic process 
with a well-defined continuous limit. 

B-S processes can be fully characterized if one regards their finite dimensional marginals 
as instances of multivariate normal variance mixtures U = aN, where a is an univariate 
random variable with positive values, having cumulative distribution G, and N is an 
77,-dimensional normal random variable independent from a. Leaving aside trivial affine 
changes of variables, we can assume that the covariance matrix of is the identity matrix. 
By first conditioning its evaluation on the value of a, and then computing its mean over a, 
it is immediate to see that the characteristic function g^i^ki, /c2, . . . , kn) of U is 

^^(/ci, k2,..., kn) = Q(A;2 + A;2 + ■ ■ ■ + kl)^ , 

where 00-2 (s) is the Laplace transform associated to G 

poo 

(j)„2{s) = / dxe~''''dG{x). 
Jo 

As this construction is independent from n, an admissible choice for is 0(s) = 0o.2(|), where 
0(j2 is the Laplace transform associated to any random variable with positive values. 

The crucial point is that by Schoenberg's theorem in Schoenberg [46] (see also the self- 
contained discussion about normal variance mixtures in Bingham and Kiesel [9]) this family 
exhausts all the possible choices, i.e. (f){kl + A;| + ■ ■ ■ + k"^) is a characteristic function of a 
probability distribution for any n if and only if 0(s) is the Laplace transform a univariate 
random variable with positive values. 

Hence a multivariate distribution for the returns can be built in the B-S framework if 
and only if it admits a representation as a normal variance mixture. 

In passing we note that the choice of B-S in their original papers for the distribution (5) 
is indeed admissible, as in Sokolov et al. [47] it is shown that 

"^'^'^ = (ttcS^) 
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is completely monotone, hence a Laplace transform by the virtue of Bernstein's theorem. 

Now it is immediate to see that all the stochastic processes X"{(jj) that can arise in the B-S 
framework admit the following representation on a suitably chosen stochastic basis (^2, JF, P), 
over which a positive random variable o-{uj) and a Wiener process Wtioj) independent from 
a are defined: 

X^{uj) = (T{uj)Wt2D{uj). (9) 

We only have to show that the finite dimensional marginal laws of X"{uj) are the same as 
those arising from (8). Indeed if we first evaluate the expectations over W ^ conditional on 
cr, we will obtain a Gaussian multivariate distribution 



1 



(27ra2 



exp 



t2D +2D +2D +2D +2D 



\t\^ tr-tr e-Ci 

the eventual average over a will then lead to the same multivariate normal variance mixtures 
as in (8), with the appropriate covariance matrix (just note that Atj = tj+i — tj, and 
fiAU — ~ ^ti)- In particular, the processes introduced in Sec. II correspond to an 

inverse Gamma distribution of cr^ in the Student case, and a Generalized Inverse Gaussian 
distribution in the hyperbolic case. 

The stochastic differential equation obeyed by (9) is 

This equation shows that the volatility of the processes admissible in the B-S framework 
has a deterministic time dynamic, and that its source of randomness is just ascribable to its 
initial value. 

Eventually we can conclude that a stochastic process is compatible with the B-S frame- 
work if and only if it is a variance mixture of Wiener processes whose variance is distributed 
according an arbitrary positive law, with a deterministic power law time change. This ex- 
plains why using use this framework to model real price returns, one inevitably has to assume 
that the real price dynamics is composed by sequences of different realizations, as done by 
B-S. This is necessary not only because otherwise the model would predict a persistent and 
deterministic volatility decay for D < 1/2, but also because a is fixed in each realization. 
The limitations of this kind of models in describing real returns will be made more manifest 
in the following section, but now we already know their mathematical foundations. 
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The asset prices can be modeled in an obvious arbitrage free way 
S{t,u) = 5*0 exp (rt + a{u;)Wt2D{u) - ^(T^{uj)t'^^ 



with r the fixed default free interest rate, and where we left the dependence on u explicit in 
order to emphasise the fact that a is a random variable. The pricing of options is then the 
same as in the Black-Scholes model, with an additional average over criuj). For instance the 
price C(T, K) of a call option with maturity T and strike K is 

C{T,K) = S,E,{N{d{)) - e-'-^KE^iNid^)) , 

with as usual N is the normal cumulative distribution. 

In ^ + rt + ^aH^^ 



K 

In f + rt - 



d. --^ ' 2 



^2 



atD 

and the additional expectation has to be evaluated according to the distribution of a. 



IV. APPLICABILITY OF THIS FRAMEWORK TO REAL MARKETS 

The axiomatic nature of the derivation of Baldovin and Stella is elegant and powerful: 
its ability to build mathematically multivariate price return distributions from a univariate 
distribution using only a few reasonable assumptions is impressive. Nevertheless, as stated 
in the introduction, a model of price dynamics must meet many requirements in order to be 
both relevant and useful. In this section, we examine its dynamics thoroughly. 



A. Volatility dynamics 

In Fig. 3. a we report the results of three simulations of the return process, each one of 
500 steps and with parameters u = 3.2 and D = 0.20. In each run the volatility decays 
ineluctably, as explained in the previous section. Indeed by fixing the time interval 6ti = 1, 
we see from Eq. (8) that the unconditional volatility of the r^ i returns is proportional to 
A/(t + - t"^^, i.e., to t^"^/^ for t > 1: the unconditional volatility decreases if D < 1/2 
and increases if D > 1/2, in both cases according to a power law. This appears quite clearly 
in Fig. 3.b, where we have computed the mean volatility decay, measured as the absolute 
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(a) Three simulations, each 500 steps long. 



(b) Decay of the volatility: average over 10000 simulation, 
each 500 steps long. The dashed line represents the 
analytic prediction. 



Figure 3: Process simulation with v = 3.2, D = 0.20, and A = 0.107. 



values of the return, over 10000 process simulations. The parameters of the distributions 
have been chosen close to those representing real returns (see below). 

The conditional volatility can be easily computed: the distribution of the return r„ i 
conditioned to the previous return realizations ro,i, . . . , Vn-i^i is again a Student distribution 
with exponent ul = u + n and conditional variance 



n-l 



[(n + 1) 



1 + V ^ 

^ fz+ 1)2^ 



i=0 



From this expression it is clear that volatility spikes in a given realisation of the process 
tend to be persistent (see Fig. 3. a); this is the main reason why fluctuation patterns differ 
much from one run to an other. This can be also understood by appealing to the character- 
ization of this kind of processes we did in Sec. Ill: each single run is just a realization of a 
Wiener process, whose variance is chosen at the beginning according to an Inverse Gamma 
distribution -Rr(|, |), and that decays in time according to the deterministic law t^~5. 
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B. Decreasing volatility and restarts 

The very first model introduced by B-S has constant volatihty, which corresponds to A 
being a multiple of the identity matrix. This unfortunate feature is the main reason behind 
the introduction of weights, whose effect is akin to an algebraic stretching of the time, or, as 
put forward by B-S, to a time renormalization. This in turn causes a deterministic algebraic 
decrease of the expectation of the volatility, as explained above and depicted in Fig. 3.b; 
hence the need for restarts, each attributed to an external cause. 

Although this dynamics may seem quite peculiar, such restarts are found at market 
crashes, like the recent one of October 2008, which are followed by periods of algebraically 
decaying volatility. This leads to an analogous of the Omori law for earthquakes, as reported 
in Lillo and Mantegna [36] and Weber et al. [53]. The B-S model, by construction, is able 
to reproduce this effect in a faithfully way. In Fig. 4 the cumulative number of times the 
absolute value of the returns N(t) exceeds a given thresholds is depicted, for a single simu- 
lation of the process and three different value of the threshold. The fit with the prediction 
of the Omori law N{t) =K{t + to)" - Kf^ is evident. 

Crashes are good restart candidates: they provide clearly defined events that synchronize 
all the traders' actions. In that view, they provide an other indirect way to measure the 
distribution of timescales of traders, which are thought to be power-law distributed (Lillo 
[35]). 

Another example of algebraically decreasing volatility was recently reported by McCauley 
et al. [40] in foreign exchange markets in which trading is performed around the clock. Under- 
standably, when a given market zone (Asia, Europe, America) opens, an increase of activity 
is seen, and vice-versa. Specifically, this work fits the decrease of activity corresponding to 
the afternoon trading session in the USA with a power-law and finds an algebraic decay with 
exponent t] = 0.35; this is exactly the same behavior as the one of B-S model between two 
restarts, with D = 1 — 21] = 0.3. No explanation of why the trading activity should result in 
this specific type of decay has been put forward in our knowledge. In this case the starting 
time of the volatility decay corresponds to the maximum of activity of US markets. 
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Figure 4: Omori law for a single run of the process, with D = 0.20, z/ = 0.32. N(t) is the 

cumulative number the absolute value of the return exceeds a given thresholds. Three 
different values of the threshold / have been chosen, measured with respect to the standard 
deviation a of the data. The dashed lines represents the fit with the Omori law 



C. Apparent multifractality 

The Baldovin and Stella model is able to reproduce the apparent multifractal character- 
istics of the real returns, i.e. the shape oi ({q) where (Ir^^l'^) = 

The expectation is evaluated according the distribution (4), i.e. taking the mean over 
independent runs of the process. Hence the expectation of the gth moment in this model is 



(see the addendum to Baldovin and Stella [5]). The exponents ({q) are evaluated as the 
slopes of the linear fitting of ln((|r|'^)p^J with respect to \ia{St). Hence in our case they are 
determined by the expression ln^[Z!^^[(t + (5t)^^ — t^^]^/^, and depend only on D and Tc- In 
Fig. 5. a is depicted the fitting of the S&P 500 exponents with the model (10). The best fit 
is obtained with D = 0.212 and Tc = 5376. Unfortunately a value of Tc that large is difficult 
to justify, as in the case of S&P 500 we have only 14956 daily returns, i.e. less than three 



N{t) = Kit + to)" - Kt\ 




(10) 
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Scaling exponents 
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(a) Fitting of the empirical exponents of real data. (b) Theoretical prediction compared to 5 simulations done 

with the same parameters. 

Figure 5: Scaling exponents: S&P 500 data and simulations compared with theoretical 
prediction. All the simulations have been done with the same parameters: 30 runs of 500 

steps, with u = 3.2, D = 0.220 



runs of a process with such a length. The other fit is obtained by first fixing Tc = 500, as in 
Baldovin and Stella [5] and yields D = 0.220. 

The statistical significance of this approach seems anyway questionable. In Fig. 5.b we 
compare the theoretical expectation of the exponents with simulations. We choose the 
parameters Tc = 500, D = 0.220 both for simulations and analytic model, with u = 3.22. 
The number of restarts in the simulation is 30 in order to have a number of data points 
similar to the S&P 500. It is evident that the exponents evaluated from the simulated data 
have a really large variance. 

The problem is that if the tail exponent u = 3.22, from an analytic perspective the mo- 
ments with q > 3.22 are infinite, hence, should not be taken into account in the multifractal 
analysis (for an analytic treatment of multifractal analysis see Jaffard [32, 33], Riedi [45]). 
The situation is somehow different in the case of multifractal models of asset returns (Bacry 
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et al. [2], Mandelbrot et al. [39]), where the theoretical prediction of the tail exponents of the 
return distribution is relatively high (see the review of Borland et al. [13]), and the moments 
usually empirically measured do exist even from the analytic point of view. For attempts 
to reconcile the theoretical predictions of the multifractal models with real data see Bacry 
et al. [3] and Muzy et al. [42]. 

It is worth remembering that the anomalous scaling of the empirical return moments 
does not imply that the return series has to be described by a multifractal model, as already 
pointed out some time ago in Bouchaud [14] and Bouchaud et al. [16]: the long memory of 
the volatility is responsible at least in part for the deviation from trivial scaling. A more 
detailed analysis of real data reported in Jiang and Zhou [34] seems indeed to exclude evident 
multifractal properties of the price series. 

V. MISSING FEATURES 

Since in this model the volatility is constant in each realization and bound to decrease 
unless a restart occurs, it is quite clear that it does not contain all the richness of financial 
market price dynamics. Restarting the whole process is not entirely satisfactory, as in 
reality the increase of volatility is not always due to an external shock. Volatility does 
often gradually build up through a feedback loop that is absent from the B-S mechanism. 
Thus, large events and crashes can also have a endogenous cause, e.g. due to the influence of 
traders that base their decisions on previous prices or volatility, such as technical analysts 
or hedgers. A quantitative description of this kind of phenomena is attempted for instance 
in Sornette [48], Sornette et al. [49], by appealing to discrete scale invariance (see also the 
viewpoint expressed in Chang and Feigenbaum [19] and references therein). This kind of 
effect is completely missing from the original B-S mechanism. 

Volatility build-ups can be simulated with D > 1/2, getting at constant D the equivalent 
of the inverse Omori law for earthquakes [29]. This kind of dynamics has been reported 
to happen prior to some financial market crashes [49]. At a smaller time scale, foreign 
exchange intraday volatility patterns have a systematically increasing part whose fit to a 
possibly arbitrary power-law, as performed in McCauley et al. [40] (77 = 0.22), corresponds 
indeed to choosing D = 0.56. To our knowledge, volatility build-ups either do not follow a 
particular and systematic law, or perhaps have not yet been the objects of a thorough study. 
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Because of the symmetric nature of all the distributions derived above, all the odd mo- 
ments are zero, hence, the skewness of real prices cannot be reproduced. This shows up well 
in Fig. 3 of Baldovin and Stella [6]. Another consequence is that it is impossible to replicate 
the leverage ejfect, i.e. the negative correlation between past returns and future volatility, 
carefully analyzed in Bouchaud et al. [17]. 

In any case, the decrease of the fluctuations in the B-S process is a deterministic outcome 
of the anomalous scaling law with D < 1/2, and results in a strong temporal asymmetry of 
the corresponding time series. But quite remarkably it misses the time-reversal asymmetry 
reported in Lynch and Zunibach [37] and Zumbach [55]. Indeed real financial time series 
are not symmetric under time reversal with respect to even-order moments. For instance, 
there is no leverage effect in foreign exchange rates, and their time series are not as skewed 
as indices, but they do have a time arrow. One of the indicators proposed in Lynch and 
Zumbach [37] is the correlation between historical volatility cr^tli^) realized volatility 
(j'fj^it). The historical volatility series <yg^^^(t) represents the volatility computed using the 
data in the past interval [t — 6th, t], and cr^t^(t) represents the volatility computed using the 
data in the future interval [t,t + 6tr]; the correlation between the two series is then analyzed 
as a function of both 6tr and 6th- Real financial time series present an asymmetric graph with 
respect the change 6th 6ts, with a strong indication that historical volatility at a given 
time scale 6th is more likely correlated to realized volatility with time scale 6tr < 6th, with 
peaks of correlation at time scales related to human activities. The asymmetry characteristic 
is absent in the Baldovin and Stella model, as showed in Fig. 6. 

The strong correlation between returns guarantees the slow decay of the volatility but 
induces some side effects. The distribution of the returns in the model is essentially the 
same with identical power law exponent for the tails. This happens independently of the 
time interval 6t over which the returns are evaluated, as long as 6t <^ Tc, with Tc of the order 
of hundreds days. Hence the weekly returns are distributed as the daily returns, while in 
real data the tail exponent begins to increase in a remarkable way already at the intraday 
level (Drozdz et al. [25]). The strong correlation also slows down the convergence to the 
Gaussian distribution of the returns when measured on larger time scale. Even if the kurtosis 
is not defined analytically in principle, it is possible to measure the empirical kurtosis of the 
returns of a simulated time series and compare with the kurtosis of real data. In Fig. 7 we 
show the kurtosis of the return distribution among simulations and daily return of the S&P 
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(a) (b) 

Figure 6: Correlation between historical and realized volatility of the simulated process, 
over different time interval 6t. The analyzed time series was composed by 1000 runs of the 
basic process, each one with 200 steps, and parameter u = 3.22, D = 0.20. 

500 index; the kurtosis has been computed for the returns over different interval 6t, and the 
simulated processes had the same length (30 runs of 500 steps) of the real series. 

VI. SUGGESTED IMPROVEMENTS 

The main limitations of the model proposed by Baldovin and Stella are poor volatility 
dynamics, lack of skewness, some unwanted symmetry with respect to time, and extremely 
slow convergence to a Gaussian. In this final section we put forward briefiy some qualitative 
proposals of how these issues can be addressed. 

The volatility dynamics can be improved by introducing an appropriate dynamics for 
the exponent D, i.e. introducing a dynamic D(t) controlling the diffusive process. This 
is equivalent to starting with a model with constant volatility, i.e. with A proportional to 
the identity matrix, and then introducing an appropriate evolution for the time t. This 
technique is employed for instance in the Multifractal Random Walk model (Bacry et al. 
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Figure 7: Comparison of the kurtosis of the returns evaluated over a time interval 5t. Each 
one of the three simulations are composed by 30 runs, 500 steps long, in order to have a 
length comparable with that of the S&P 500 returns. The parameters are u = 3.2, 



[2]), where the time evolution is driven by a multifractal process, or when the time evolution 
is modeled by an increasing Levy process (see e.g. Cont and Tankov [22]). In this last case 
we would obtain a mixing of Wiener processes driven by a subordinator. 

The lack of skewness is a common problem of stochastic volatility models: one usually 
writes the return at time t as rt^st = ^(t)<^{'t)j where e{t) is sign of the return and a(t) its 
amplitude, a symmetric setting if the distribution of e(t) is even. One remedy found for 
instance in Eisler and Kertesz [27] is to bias the sign probabilities while enforcing a zero 
expectation; more precisely, 



Another possibility for introducing skewness is that of considering normal mean-variance 
mixtures, instead of simply normal variance ones. For instance, this would have implied the 
use of the multivariate skewed Student distribution in the model described in Sec. II. 

The decay of the tail exponent of the return distribution, represented in Fig. 7, could be 



D = 0.20, A = 0.1. 
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implemented by introducing two different Student distributions: a univariate witli exponent 
Ur for modeling the daily returns, and a multivariate one with a much larger exponent 
for modeling the correlations among them. By taking into account the generalized central 
limit theorem expressed in Eq. (3), the distribution of returns at intermediate time scales 
will interpolate between the two exponents, yielding the desired feature. 

The Zumbach mugshot is one of the most difficult stylized facts to reproduce. To our 
knowledge the best results in that respect was achieved in Borland and Bouchaud [12], where 
a specific realization of a quadratic GARCH model is introduced, motivated by the different 
activity levels of traders with different investment time horizons, which take into account 
the return over a large spectrum of time scales. More specifically Borland and Bouchaud 
use 



2 2 



1 + E 



2 



gAt- 



5t 



5t=l 



with r fixing the time scale, rt sT = lnS'(t + 6T) — ln5'(t), gst measuring the impact on the 
volatility by traders with time horizon 5t, and chosen by the authors gst = 9 /{St)". This 
expression is rewritten also in the form 



j<i,k<i 

with 

oo 

M{^,J,k)= Yl f- 

At=max(i— fc) 

In the present framework this would correspond to use a highly non-trivial matrix A, 
introducing linear correlation among returns at any time lag. This means that the B-S 
process would no longer be a model of returns, but of stochastic volatility. 



VII. DISCUSSION AND CONCLUSIONS 



When employed with self-decomposable distributions like the Student or the Generalized 
Hyperbolic as introduced in Sec. II, the resulting description of the process return is different 
than that of other models in the literature. First our Student process is not stationary, hence 
different from the class of Student processes discussed in Heyde and Leonenko [30], where the 
main focus is on stationary ones. The processes (9) are also different from the one studied 

19 



in Borland [11]: the latter too are continuous and based on the Student distributions, but 
defined by the stochastic differential equation 



apart from the striking difference with Eq. (9), in Vellekoop and Nieuwenhuis [51] it is shown 
that not all the marginal distribution laws of Xt are of Student type. 

Instead in Eberlein and Keller [26] the Generalized Hyperbolic laws are adopted for 
describing the returns at a fixed time scale; these laws are then extended to the other time 
scales using the standard Levy process construction: in this case the distributions at the 
other time scales are no more of Generalized Hyperbolic type. 

The Baldovin and Stella model is also intrinsically simpler than the ones described in 
Barndorff-Nielsen and Shephard [7], where the volatility has a dynamic modeled by Ornstein- 
Uhlenbeck type processes, 



driven by an arbitrary Levy process Lt. In this case, according to the choice of L^, any self- 
decomposable distribution (like the Generalized Inverse Gaussian, or any of its special cases, 
like the Inverse Gamma) can arise as the distribution of a"^ for any t. But this simplification 
comes at a high price: while in Barndorff-Nielsen a is truly dynamic, it is fixed in B-S for 
any single process realization. 

In addition, the models analyzed in Carr et al. [18] are of a different type, even if there 
are some analogies in the underlying principles. In Carr et al. [18] indeed an anomalous 
scaling is introduced by considering self-similar processes, and in that framework any self- 
decomposable distribution can employed for modeling returns, but once again only at a 
fixed time scale, as in the standard case of Levy processes. The main difference is that in 
Carr et al. [18] the returns at different times are assumed to be totally independent, but 
not identically distributed: instead Baldovin and Stella assume that the returns are only 
linearly independent, but now with identical distributions at all the time scales, up to a 
simple rescaling. 

In conclusion, despite its current inability to reproduce all the needed stylized facts, the 
new framework proposed by Baldovin and Stella introduces a new mechanism for modeling 
returns, based on a few reasonable first principles. We therefore think that, once suitably 




t 



Xa^dt + dLt 
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modified for instance along the lines proposed above, the B-S framework can provide a new 
tool for building models of financial price dynamics from reasonable assumptions. 



Appendix: Some Useful Facts About Student and Symmetric Generalized Hyper- 
bolic Distributions 



Characteristic function of Student distributions 



The standard form of univariate Student distribution is 

r(| + i) 1 



while the multivariate one is 



T(- + - 



W2r(|) (1 +r2)f+i 



with r = ^/Y.7=l A and V{r > R) (x l/W. 

Using some standard relationships involving Bessel functions one can compute analyti- 
cally the corresponding characteristic function: 



/ + CXD 
dxi e^^i^i^i(xi) 
-oo 

'^KV Jo r(2) 

with k = |/ci|, Ka the modified Bessel function of third kind, and the employ of identity 
7.12.(27) of Erdelyi [28] 

K,{z) = ^7^^(z/ + -) / dt {t^ + z^-"-'/^ cos(t) 

1 TT 

3fJ(z/) > --, I arg(2;) |< -. 

For an alternative derivation we refer to Hurst [31] and to the discussion in Heyde and 
Leonenko [30]. An alternative expression is found in Dreier and Kotz [24]. 
For general n we obtain again the same expression. Indeed 
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r(f + f} ''+^ 
W2r(f; 

on/2'nf u+n \ r+oo 

= ^ k'-'^ / rfrr"/2(l + r^)-S-5j„,/2-i(fcr) 



V n 
2 



with = \/Yl^=i ki, c?"^^r2 the surface element of the sphere 5*"^^, the angle between k 
and X and the employ of identities 7.12.(9) 



^i^)>-l, (11) 



and 7.14.(51) of Erdelyi [28], 



^{2u -h> 3?(/i) > -1, ^{z) > 0. 

Eventually one finds 

^„(k) = ^1 (^y^/cf + --- + A;2^ . 

With the linear change of variables x C^^x, setting = (C-^)~^C~^, i.e. A = CC^, 
one obtains the following generalizations: 

rf- + -) 1 

9n{^) W2(det A)i/2r(f ) (1 + x*A-ix)f ' ^ ' 

with characteristic function 

^fn(k) = |;^(k*Ak)ti^.((k*Ak)V2). 

In the univariate case A is substituted by the scalar and the previous expressions 
reduce to 

^i(^) = li/KwIw. (13) 



'£ _|_ 1 

vrV2Ar(|) (l + f|)f+i 
and 



gi{k) = '^{Xk)^K.{\k). 
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Moments of Student distributions 



Due to the symmetry under reflection all the odd moments vanish. For the second 
moments we have, provided that z/ > 2, 

A,, 



u-2 

The moments of order 2n exist provided that z/ > 2n ; as happens for Gaussian distributions, 

they can be expressed in term of the second moments, 

r(f-n) -|-|. 

2 ' all the pairings 

In the univariate case these formulas reduce to Eix'^) = and 

(2n- l)!!r(^ -n) ^ 

' 2«r(f) ^ ■ 

u-2 



The kurtosis is then k = 3^3|, provided that u > 4. 



Simulation of multivariate Student distributions 

The simulation is a standard application of the technique used in the case of rotational 
invariance. From 

with r > 0, we see that the density of the angular variables is uniform, while setting y 



l+r2 ' 



with I > y > and r = y^y/{l — y), the density of y is given by 

1 



^l2' 2 



-?/2 (l-y)2 dy, 



i.e. by the beta distribution with parameters | and |. Eventually we can simulate the 
multivariate n dimensional distribution by 



1. Simulating y according to |) and setting r = y yz^- 

2. Simulating n i.i.d. Gaussian variables Ui and settings n 



3. Returning xn. 

The more general case (12) is simulated using the same algorithm and then returning Cx, 
where A-^ = {C^Y^C~^ , i.e. A = CC^. 
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Characteristic function of symmetric generalized hyperbolic distributions 



We start from the expression 



f (^\ ^ ^ 2 + 2^ I 

^"^ ^ (27r)tir.(a) (l + r2)f+f ' 

with r = a/XIj-Li the general case is obtained simply with an affine transformation 
X /i + 5Rx, with |U G M", 5 > a scale parameter, and R an orthogonal transformation in 
W . The central expression we need is an integral of the Sonine-Gegenbauer type, cf. identity 
7.14.(46) of Erdelyi [28]: 



/■OO 

Jo 



3fJ(/i) > -1, ^{z) > 0. 

For n = 1, considering that J_i(x) = \/ — cos{x), we obtain 



fi{ki) = I rfxie'*^^"i/i(xi) = ^ ^ / cixi \ 1 ^' cos(A;iXi) 

(27r)2K.(a) Jo (l + xf)4 + 4 



Jo 



dxiJ_i{kiXi) K^j^i{aJl + xl){l + Xi)~^~^a;f 



For alternative derivations in the univariate case see Hurst [31] and the references therein. 

In our setting the computation is exactly the same for general n, with k = \JYll=i ^"ii 
dP'~'^VL the surface element of the sphere S*""^, (p the angle between k and x, using identity 

(11) 

/„(k)= / d"xe*^-Vn(x) 

-/"-2o / ^^^n-l / ^\^n-2( I \ikr cos (t> 2^2 ^ I £ 



a 2 



/n+00 nn 
d''-^Q drr""-^ c/0 sin"-2(0)e" 
Jo Jo 



(27r)tK|(a)J Jo Jo (l + r2)t+f 

k 2 Q/2 



KE(a) 



p+oo 

/ dr J!i_i{kr)K!^_^n{aVl + r2)(l + r^)~^"5rt 
Jo ^ ^ ^ 
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Hence the eventual result /n(k) = fi{k). 
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